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In a previous paper [J. Chem. Phys. 129, 024506 (2008)] we studied a 3-dimeiisionaI lattice 
model of a network-forming fluid, recently proposed in order to investigate water anomalies. Our 
semi-analytical calculation, based on a cluster-variation technique, turned out to reproduce almost 
quantitatively several Monte Carlo results and allowed us to clarify the structure of the phase dia- 
gram, including different kinds of orientationally ordered phases. Here, we extend the calculation 
to different parameter values and to other similar models, known in the literature. We observe that 
analogous ordered phases occur in all these models. Moreover, we show that certain "waterlike" 
thermodynamic anomalies, claimed by previous studies, are indeed artifacts of a homogeneity as- 
sumption made in the analytical treatment. We argue that such a difficulty is common to a whole 
class of lattice models for water, and suggest a possible way to overcome the problem. 



I. INTRODUCTION 

Water is extremely abundant in nature, and it is of 
enormous importance from a biological, technological as 
well as environmental point of view, because of its var- 
ious anomalous properties For instance, it is well 
known that liquid water exhibits unusually large heat 
capacity and dielectric constant. Furthermore, at ordi- 
nary pressures, the solid phase (ice) is less dense than 
the corresponding liquid phase, while the latter displays 
a temperature of maximum density, slightly above the 
freezing transition. In spite of great research efforts on 
water [1, S I3|> ^ f^l^y consistent theory of such (and 
many other) anomalies from first principles is not yet 
available. Nonetheless, it is well established that most 
anomalies are related to the ability of water molecules to 
form a network of hydrogen bonds. 

Following this view, so-called network-forming fluids 
have been the subject of several theoretical investiga- 
tions. Among different possible approaches, a number 
of studies have been developed m the framework of sim- 
plified lattice models 
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Various types of model molecules with 



orientation-dependent interactions have been used, in ei- 
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mensions. A natural choice for water is a 3-dimensional 
model molecule with four bonding arms arr ang ed in a 
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rietrpiilJi^ UJ, [14|, [151, ll6|, ll7|, ll8|, ll9|, 

120,, ll, 12, 23, 2A, 113 . Two arms represent the hy- 

drogen (H) atoms, which are positively charged and act 
as donors for the H bond, whereas the other two repre- 
sent the negatively charged regions of the H2O molecule 
("lone pairs"), acting as H-bond acceptors. As far as the 
lattice is concerned, the body-centered cubic (bcc) lattice 



is suitable for the tetrahedral molecule, as the latter can 
point its arms toward four out of eight nearest neighbors 
of each given site. The above features are common to 
several models, which differ in the form of interactions 
and in the set of allowed configurations. 

In the early model proposed by Bell [ill, [13, [ll], 
molecules can point their arms only toward nearest neigh- 
bor sites. An attractive energy is assigned to every pair of 
occupied nearest neighbors, with an extra contribution if 
a H bond is formed, i.e., if a donor arm points toward an 
acceptor arm. Moreover, a repulsive energy is assigned to 
certain triplets of occupied sites, in order to account for 
the difficulty of forming H bonds by closely-packed water 
molecules. Minor variations of this model have also been 
investigated: Bell and Salt [l3|, and subsequently Meijer 
and coworkers [l^ [l^ , have replaced the three-laody in- 
teraction with a simple next-nearest-neighbor repulsion, 
whereas Lavis and Southern [l7| have defined a simpli- 
fied model with no distinction between donors and accep- 
tors. All these studies predict a liquid- vapor coexistence 
and two different low-temperature phases, characterized 
by orientational order and different densities. At zero 
temperature, the low-density phase is an ideal diamond 
network made up of H-bonded water molecules, with half 
the lattice sites left empty, resembling the structure of ice 
Ic (cubic ice) . Conversely, the high-density phase is made 
up of two interpenetrating diamond structures, with all 
the sites occupied, resembling the structure of ice VII 
(a high pressure form of ice) . If the models retain bond 
asymmetry, i.e., donors and acceptors are distinguished, 
then both zero-temperature phases possess a residual en- 
tropy. 

Debenedetti and coworkers have studied a similar 
model ^Sk [m, [111 , in which water molecules have an ex- 
tra number of nonbonding configurations and the close- 
packing energy penalty occurs only when two molecules 
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in a triplet form a H bond. At zero temperature, this 
model still predicts the two ordered H-bond networks of 
the Bell model. At finite temperature, these (icelike) 
phases have not been investigated, as the cited works 
were mainly focused on metastable liquid water. Indeed, 
the model seems to support the popular "second critical 
point" conjecture ]33, |33|, originally proposed by Stan- 
ley and coworkers [34| . Two authors of the present paper 
have also studied a simplified version of this model, with- 
out the donor-acceptor asymmetry [H, [2I] . 

A variation of the Bell model has been considered by 
Besseling and Lyklema [l^, [l^, who have taken into 
account only nearest-neighbor interactions, namely, an 
attractive term for H-bonded molecules and a repulsive 
term for nonbonded molecules. Even in this case, the 
two different ordered phases described above are stable 
at zero temperature, but they have not been investigated 
at finite temperature. The cited studies were indeed de- 
voted to liquid-vapor interface properties and to hy- 
drophobic hydration thermodynamics [l^ , for which the 
authors found good agreement with experiments. These 
results were obtained in the so-called quasi-chemical or 
Bethe approximation, i.e., a first-order approximation 
which takes into account correlations over clusters made 
up of two nearest-neighbor sites. 

Girardi and coworkers [2^ have recently performed ex- 
tensive Monte Carlo simulations for the above model, in 
the simplified version with symmetric bonds. This work 
claims the onset of two liquid phases of different densities, 
which can coexist in equilibrium, in agreement with Stan- 
ley's conjecture [s^l- The coexistence line seems to ter- 
minate in a critical point, while the lower-density phase 
exhibits a temperature of maximum density, depending 
on pressure. 

In our previous paper [23| we analyzed the latter model 
by a generalized first-order approximation, based on a 
four-site tetrahedral cluster. Such semi-analytical calcu- 
lation turns out to reproduce, with remarkable accuracy, 
most physical properties obtained by simulations. Nev- 
ertheless, the resulting phase diagram is quite different 
from the one proposed in the original paper [S^. The 
two different condensed phases exhibit orientational or- 
der, which suggests to identify them as a temperature 
evolution of the two zero-temperature network structures 
introduced above. Furthermore, the two claimed critical 
points turn out to be in fact tricritical, and two critical 
lines appear, related to two different kinds of symmetry 
breaking. In conclusion, the phase diagram is more com- 
plex and richer than expected, but unfortunately quite 
far from the real water phase diagram. 

In the current paper, we exploit the effectiveness of our 
approximation scheme, with a twofold purpose. On the 
one hand, we analyse in more detail the simple model by 
Girardi and coworkers (GBHB), exploring the phase dia- 
gram for different parameter values. For strong repulsive 
interaction (slightly weaker than the H bond) , we obtain 
an even richer phase diagram, including a new ordered 
phase, appearing at finite temperature. On the other 



hand, we reconsider more complex models with asym- 
metric bonds, such as the Besseling-Lyklema (BL) model 
and the original Bell model. 

Several investigations devoted to this kind of mod- 
els have hypothesized a fully homogeneous (i.e., site- 
independent) probability distribution of molecular con- 
figurations, thus excluding the possibility of orientational 
ordering of the type observed in the GBHB model [l^ . 
Conversely, we find out that analogous ordered phases 
occur even in the BL and Bell models. The liquid phase 
reported by previous papers [ill, [3 (displaying water- 
like density anomalies) turns out to be an artifact of the 
aforementioned homogeneity assumption. Furthermore, 
density anomalies appear in the lower-density ordered 
phase, which makes it questionable to regard this phase 
as a representation of some real ice form. Let us note that 
density anomalies in an icelike phase have been previ- 
ously noticed by Meijer and coworkers [T6| , investigating 
the Bell-Salt model [14|. We discuss the reasons under- 
lying these difficulties, which we argue to be common to 
a wide class of network-forming lattice models for wa- 
ter, and finally suggest a possible way to circumvent the 
problem. 

The paper is organized as follows. In Section II we 
introduce a generic model hamiltonian incorporating the 
three models under investigation. In Section III we de- 
scribe in more detail the ground states of the models. 
In Section IV we explain the cluster-variation technique 
employed for the calculation. Section V reports the re- 
sults, and a comparison to previous semi-analytical [Tsj 
and numerical [13] investigations. Section VI contains 
the conclusions. 



II. THE MODEL HAMILTONIAN 

As mentioned in the Introduction, we study three dif- 
ferent models of increasing complexity (GBHB, BL, and 
Bell), defined on a bcc lattice (Fig. [T]). It is possible to 
write a single hamiltonian incorporating the three mod- 
els. For all the models, a lattice site can be empty or 
occupied by a molecule having a tetrahedral structure (4 
bonding arms, which can point toward 4 out of 8 nearest 
neighbors of each given site). 

In the simplest case (GBHB model), a hydrogen bond 
is formed, yielding an attractive energy —77 < 0, when- 
ever two nearest-neighbor molecules point an arm toward 
each other, with no distinction between donor and ac- 
ceptor. Each molecule can thus have only two different 
configurations, which we simply denote as 1 and 2 (see 
Fig. [2]). Moreover, a repulsive energy — e > is assigned 
to every pair of nearest-neighbor sites occupied by wa- 
ter molecules. This term is meant to penalize neighbor 
molecules not forming H bonds. Finally, a chemical po- 
tential contribution —fj, is taken into account for every 
occupied site (in this paper we always consider a grand- 
canonical description) . 

In the BL model [18], the above picture is a little more 
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FIG. 1: Two conventional (cubic) cells of the body centered 
cubic (bcc) lattice. A,B,C,D denote four interpenetrating 
face-centered cubic (fee) sublattices. 




FIG. 2: Two model molecules in a bonding configuration. 
The arm configuration is i = 1 for the lower molecule and 
i = 2 for the upper one (see the text). 



complex. The bonding arms of each molecule are "deco- 
rated" by a sign (two and two "— "), respectively de- 
noting donors and acceptors. A H bond is formed only if 
the two bonding arms, pointing toward each other, carry 
one sign and one "— " sign, or vice versa. As a con- 
sequence, a water molecule can assume 6 distinguishable 
"sign configurations" for each of the 2 "arm configura- 
tions" described above, adding up to 12 different config- 
urations. 

The Bell model is equivalent to the BL model 



in the bonding mechanism, but incorporates different 
orientation-independent interactions. The two-molecule 
interaction is attractive (— e < 0), whereas a repulsive 
energy —7 > is assigned to any (fully occupied) 3-site 
cluster, made up of two second neighbors plus one of their 
common first neighbors. 

We can write the general model hamiltonian as 

H = -e ^ rii^nj^ - rii^ (1) 

{r,s) r 
(r,s) (r,s,i) 

where r,s,t denote lattice sites, irjisjit denote their re- 
spective configurations, rii is an "occupation function", 
defined as = if i = (empty site), rii — 1 other- 
wise (occupied site) , and hij is a "bond function" , de- 
fined as hij = 1 if the pair configuration repre- 
sents a H bond, and hij = otherwise. The summations 
denoted by r, (r, s), and (r, s,t) respectively run over 
sites, nearest neighbor pairs, and the 3-site clusters de- 
fined above. Let us note that the first line of Eq. ([T]) 
is formally equivalent to a lattice-gas hamiltonian, al- 
though the GBHB and the BL models are characterized 
by a repulsive interaction energy (e < 0). The first term 
of the second line represents the H-bond energy, which 
is always attractive (77 > 0), whereas the last term takes 
into account the 3-body repulsive interaction (7 < 0). As 
mentioned above, the latter occurs only in the Bell model, 
whereas the GBHB and BL models are characterized by 
7 = 0. Concerning configuration indices, we can always 
denote a vacancy (empty site) by z = 0, but molecu- 
lar configurations need to be defined in different ways, 
depending on whether the model assumes asymmetric 
bonds (donor-acceptor distinction) or symmetric bonds 
(no distinction). In the latter case (GBHB model), the 
two possible arm configurations, which we respectively 
denote by i = 1, 2 (see Fig. [2]), provide a full description 
of the molecule configurations. In the other case we have 
to understand that i is in fact a multi-index (i, i'), where 
i' denotes the sign configuration. As a consequence, in 
general the precise definition of the bond function hij 
turns out to depend on the spatial orientation of the lat- 
tice bond considered. It can be made unambiguous in the 
symmetric bond case, by defining a special order for the 
sites in each nearest neighbor pair. Such order may be 
for instance the one indicated by the arms of a molecule 
in the i — 1 configuration, in which case we obtain the 
following simple definition: hij = 1 if i = 1 and j = 2, 
and hii = otherwise. 



III. GROUND STATES 

The models under investigation can exhibit two dif- 
ferent ground states with different densities, in which 
the bonding arms of the molecules form ordered net- 
work structures, as described in the Introduction. Such 
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ordered structures naturally split the bee lattice into 
four (mutually exclusive) interpenetrating fcc sublattices, 
which we denote by A,B,C,D (see Fig. [1]). Assuming 
here that i, j,k,l denote the arm configurations of all sites 
placed on the A, B, C, D sublattices, respectively, one can 
see that the low density structure (single diamond net- 
work) is fourfold degenerate and can be represented by 
the four alternative sublattice configurations (i,j,k,l) = 
(1, 2, 0, 0), (0, 1, 2, 0), (0, 0, 1, 2), (2, 0, 0, 1) (which can be 
obtained from one another by a circular permutation). 
In each structure, two sublattices (respectively, AB, BC, 
CD, and DA) are occupied, while the other two sub- 
lattices are empty. Conversely, the high-density struc- 
ture (two intertwined diamond networks) turns out to 
be twofold degenerate and can be represented by the 
two alternative sublattice configurations (i,j, fc, Z) = 
(1, 2, 1, 2), (2, 1, 2, 1). Both these configurations have all 
lattice sites occupied, but they differ in the pairs of 
bonded sublattices, which is, AB and CD, in the for- 
mer case, or alternatively BC and DA in the latter. 

In this paper, we do not investigate in detail the ground 
state phase diagrams, which depend on the interaction 
parameters of each different model. Let us only re- 
mark an important difference between the simpler GBHB 
model, with symmetric bonds, and the other two mod- 
els with asymmetric bonds. In the former case, molecule 
configurations are fully defined by their arm configura- 
tions, so that the formation of a diamond network imme- 
diately implies that all molecules are maximally bonded 
(each molecule can participate in four bonds at most). 
Conversely, in the latter case, molecules still have the 
freedom of arranging their sign configurations in order 
to respect the ice rule (i.e., to realize the correct donor- 
acceptor pairings), giving rise to an exponential number 
of global configurations with equal energy, and there- 
fore to a residual ground-state entropy. An approxi- 
mate calculation of this entropy was first performed in 
1935 by Pauling, who obtained a value of ln(3/2) per 
molecule [S^. This value turns out to be in extremely 
good agreement with both experimental measurements 
and simulations (see Ref. [s^ and references therein). 



IV. CLUSTER- VARIATION APPROXIMATION 

The ground state structures described above suggest 
to write the average (grand-canonical) energy per site in 
the following form 



w = Pijki'Hijku 

i,j,k,l 



(2) 



where Pijki denotes the joint probability distribution for 
the configurations of a tetrahedral cluster like the one 
depicted in Fig. [3l whose sites are placed respectively on 
the four different A, B, C, D sublattices, and TCijki is a 
suitable function, which we denote as tetrahedron hamil- 
tonian. The sum is understood to run over all possible 




FIG. 3: Tetrahedral cluster, made up of four sites placed on 
the four different A,B,C,D sublattices. AB, BC, CD, and 
DA are nearest- neighbor pairs; AC and BD are next-nearest- 
neighbor pairs. 



configurations of the four sites. By grand-canonical en- 
ergy, we mean w = u — fip, where u is the internal energy 
per site and p is the density, i.e., the average occupation 
probability 



i,j,k,l 



Hi + rij -\- Tik -\- ni 



(3) 



At first, let us consider the model with symmetric bonds 
(GBHB model) , so that the configuration indices denote 
only arm configurations. It is easy to verify that the 
tetrahedron hamiltonian can be written as follows, 



Ti-ijki = Ti-ijk + Ti-jki + Ti-kii + Huj, 



where 



(4) 



(5) 



Eq. (Ill) shows that the tetrahedron hamiltonian turns out 
to be invariant under circular permutation of the config- 
uration indices i,j,k,l. As far as Eq. ([5]) is concerned, 
numerical coefficients are meant to adjust counting of the 
different terms. The one-site (chemical potential) term is 
divided by 4, which is, the number of sites in the tetrahe- 
dron. The two-site terms have a unit coefficient, because 
there are exactly 4 nearest-neighbor pairs in the tetrahe- 
dron and 4 nearest-neighbor pairs per site in the lattice. 
Finally, the three-site term is multiplied by 3, because 
there are 4 triangle clusters in the tetrahedron, but 12 
triangle clusters per site in the lattice. 

Let us note that the meaning of Eq. and of the 
related definition of the tetrahedron hamiltonian, is that 
every elementary tetrahedron is assumed to give on aver- 
age the same contribution to the energy of the system. In 
other words, we assume that possible breaking of transla- 
tional invariance in the thermodynamic state of the sys- 
tem can only occur at the level of the four sublattices 
identified by the ground states, whereas the thermody- 
namic state of each individual sublattice is assumed to 
remain translationally invariant. This assumption is sup- 
ported by Monte Carlo simulations, performed for the 
GBHB model [13 . 

As far as entropy is concerned, we make use of 
the very same approximation employed for the GBHB 
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model [27|. In the cited paper, we deduced the ap- 
proximation as an instance of Kikuchi's cluster-variation 
method [s^, [H, [s^. The latter is a generalized mean- 
field theory, which describes correlations up to the size of 
certain maximal clusters. In general, one obtains a free- 
energy functional in the cluster probability distributions, 
to be minimized, according to the variational principle 
of statistical mechanics. Here, we derive the variational 
free energy, according to a heuristic argument analogous 
to the one originally suggested by Guggenheim [i^l • One 
assumes that the entropy per site can be evaluated as a 
difference of two terms. The former is the information 
entropy associated with the 4-site probability distribu- 
tion pijki, which takes into account correlations among 4 
configuration variables on a tetrahedral cluster (i.e., on 
the 4 different sublattices). The latter can be viewed as 
a correction term, which ensures that, if the tetrahedron 
distribution factorizcs into a product of single-site proba- 
bilities, the mean- field (Bragg- Williams) entropy approx- 
imation is recovered. 

The variational grand-canonical free energy per site 
Lu = w — Ts (T being the temperature, expressed in en- 
ergy units, and s the entropy per site, in natural units) 
can be finally written as 



Pijkl 



n 



■ijkl 



T 



(6) 

where pf^ is the probability of the i configuration for a site 
on the X sublattice {X = A, B,C, D). These probabili- 
ties are of course defined as marginals of the tetrahedron 
distribution as 



pf 



Pi 



k,l,i 

'^Ptjki- 



(7) 



The variational free energy in Eq. ([6]) is thus a function 
of the only tetrahedron distribution pijki- Let us note 
that such a free energy is also sometimes referred to as 
(generalized) first-order approximation (on the tetrahe- 
dron cluster), and turns out to be exact for a suitable 
Husimi lattice (4ll [4^ . with the (tetrahedral) building 
blocks arranged as in Fig. 31 

The free energy minimization, taking into account the 
normalization constraint 



i,j.k.l 



1, 



(8) 



can be performed by the Lagrange multiplier method, 
yielding 



Pijki — t, e [p^ Pj PkPi j 



(9) 



where C is related to the Lagrange multiplier, and can 
be determined by imposing the constraint Eq. ([5]). One 




FIG. 4: Local structure of a Husimi lattice, made up of 
tetrahedral blocks. 



obtains 



c = 



i,j,k,l 



^-/^ {pfpfp'^pFf 



(10) 



Eq. dH), together with Eqs. ^ and 0, provides 
a fixed-point equation for the tetrahedron distribution 
Pijkl, which can be solved numerically by simple iter- 
ation (natural iteration method 43]). This numerical 
procedure can be proved to lower the free energy at each 
iteration f4^.l43j. 

From the tetrahedron distribution, one can determine 
the thermal average of every observable. The density can 
be computed by Eq. ^ , the internal energy u = w + ^ip 
by Eqs. ^ and ([3]), and the free energy by Eq. ([6]). 
The latter can also be related to the normalization con- 
stant as cj = — TlnC [i^l, so that the entropy reads s = 
InC + w/T. Finally, assuming the volume per site equal 
to unity, pressure can be determined, in energy units, 
as P = —to. In the presence of multiple solutions, i.e., 
of competing phases, the thermodynamically stable one 
is selected by the lowest free energy (highest pressure) 
value. 

As far as first order transitions are concerned, they can 
be easily determined by finding (numerically) a change 
of sign in the difference between the free energy values 
of two minima of the variational free energy. Second or- 
der (critical) transitions are numerically more delicate, 
since they are characterized by a free energy minimum 
becoming a saddle point. In this case, it is convenient to 
determine changes of sign in some eigenvalue of the Hes- 
sian matrix of the variational free energy. The elements 
of the latter can be written as 



1 d^u; 
T dpijkidpi'j'k'v 
_ Sii'SjjiSkk'Sii' 3 
Pijkl 4 



(11) 



6ii' Sjji 



Skk' , Suf 



P? 



P, 



Pi 



D 



where 5 denotes Kronecker delta. 

As mentioned above, these calculations are valid in the 
case of symmetric bonds (GBHB model), but they have 
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to be generalized, in order to deal with the asymmetric- 
bond models (BL and Bell models). In the latter cases 
one has to take into account both arm and sign config- 
urations. We make the simplifying assumption that all 
sign configurations are equally probable, for a molecule 
in a given arm configuration. Such assumption is reason- 
able, since it turns out to lead to a very good approxi- 
mation for the zero-point entropy of perfect ice (where 
correlations among sign configurations are expected to 
be maximal). Indeed, in Ref. 18, Besseling and Lyklema 
showed that Pauling's calculation is exactly equivalent to 
the quasi-chemical (Bethc) approximation, with the ex- 
tra assumption of equally probable sign configurations. A 
simple, though quite tedious, calculation proves that the 
same result is recovered by our cluster approximation. In 
a few words, this is related to the fact that each tetra- 
hedral (4-site) cluster contains only 2 molecules in the 
same diamond network. Therefore, for perfect ice (i.e., 
fixed arm configuration), the tetrahedral cluster approxi- 
mation describes sign-configuration statistics in the very 
same way as the Bethe (2-site) approximation does. 

According to these arguments, we restrict the mini- 
mization procedure to a subspace of the probability dis- 
tribution space, characterized by equally probable sign 
configurations. It is possible to show that such a calcula- 
tion leads to a set of equations for the probability distri- 
bution of arm configurations, that are formally equivalent 
to Eqs. ([9]) and pO|) . In these equations, the tetrahedron 
hamiltonian Tiiju can still be evaluated by Eqs. ^ and 
but the H bond energy rj is replaced by the effective 
temperature-dependent parameter 



fy(r) Tin 



1 + e"/^ 



while the chemical potential /i is replaced by 
fi{T) ^ n + T\n6. 



(12) 



(13) 



(recall that we have 6 sign configurations for each arm 
configuration). These parameters contain entropic con- 
tributions, which take into account the extra degrees of 
freedom related to the sign configurations. Such contri- 
butions must not be taken into account in the evaluation 
of the (grand-canonical) average energy. It turns out that 
the latter can still be computed by Eqs. Q, (jH), and ([5]), 
with 7] replaced by 



r]{T) = 1] 



1 + e"/^ ■ 



(14) 



The latter parameter represents the H-bond energy, av- 
eraged over the sign configurations. 

Let us finally recall that we are also interested in study- 
ing the effect of a homogeneity constraint, i.e., of forcing 
site probability distributions pf to be independent of the 
sublattice X. It is easy to see that, if the latter condition 
is verified in the right-hand side of Eq. © the result- 
ing tetrahedron distribution (left-hand side) turns out to 
be invariant under circular permutation of the indices. 



due to the same invariance property of the tetrahedron 
hamiltonian Tiiju ■ As a consequence, the marginals com- 
puted by Eqs. ^ will verify the homogeneity condition 
as well. This observation clearly shows that the cluster- 
variational free energy actually admits "homogeneous" 
stationary points. In order to determine such points, 
preventing numerical instabilities which might drive the 
iterative procedure toward more stable symmetry-broken 
solutions, it is sufficient to compute marginals according 
to the following prescription 

A B C D \ " Pijkl + PUjk + Pklij + Pikli 

Pi =Pi =Pi =Pi =2^— — 



3,k,l 



rather than Eqs. ([7]). 



(15) 



V. RESULTS 



Phases 



Let us consider the probability distributions of the 
arm configuration i on the X sublattice, as defined in 
Eqs. ([7]). Due to normalization, each distribution can 
be represented by two parameters, which can be conve- 
niently defined as 



P^. 



(16) 
(17) 



The former is an occupation probability, restricted to 
the X sublattice, so that we may call it sublattice den- 
sity. The latter characterizes preference for the water arm 
configuration 1 rather than 2 on the X sublattice. For 
the benefit of readers that are familiar with spin models, 
let us note that the arm configuration variables of the 
current models admit a "natural" mapping onto spin-1 
variables, namely, S' = ±1 spin values representing the 
two water configurations, and S — Q representing empty 
sites. In this view, the above parameters turn out to 
be the usual order parameters of spin-1 models, i.e., the 
quadrupolar order parameter (5^) and the magnetiza- 
tion (S"), respectively. It is possible to characterize the 
different phases of our models by the set of sublattice 
order parameters p'^,^'^, for X = A, B,C, D. Let us 
stress the fact that the same characterization is valid for 
all the models under investigation, both with symmetric 
and asymmetric bonds. 

At high temperature and low pressure, we always ob- 
serve a fully homogeneous (gaslike or liquidlike) phase, 
which we simply denote as fluid (F), where the order pa- 
rameters p^,^^ are independent of the sublattice. The 
parameters vanish, revealing that there is no prefer- 
ence between the two alternative arm configurations of a 
water molecule. In summary, we can write 



p^ p° 

e e e e 



Pf Pf Pf Pf 




(18) 
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where the matrix notation has been chosen for reasons 
of clarity. In the low-temperature region, we observe 
two different phases, which we denote as high density 
(HD) and low-density (LD) phases. These phases can 
be regarded as a temperature evolution of the two (HD 
and LD) ground-state ordered network structures, re- 
spectively. 

As far as the HD phase is concerned, we have 



1.5 ■ 



HD 



(a) 



pA pB pC pD- 
e e e 



phd phd phd phd 
Chd — Chd '^hd —Cud 



, (19) 



with phd = Chd = 1 at zero temperature. The den- 
sity is still homogeneous, whereas the preference param- 
eters reveal that H bonds are more likely formed within 
the sublattice pairs AB and CD, respectively. In other 
words, this phase can be viewed as an ideal HD struc- 
ture (two interpenetrating diamond networks), in which 
defects are progressively introduced by thermal fluctua- 
tions. Like in the ground state, this phase turns out to 
be twofold degenerate. The degenerate solution, having 
the same free energy, can be obtained from Eq. p9p by 
a circular permutation of the sublattice indices. 
The LD phase is characterized by 



p^ p^ p^ p^' 



PhD PlB PhD PhD 
^LD ^?LD ClD ~ClD 



(20) 



Assuming for instance p^u > Pld i the sublattice pair AB 
turns out to be more populated than CD, whereas the 
preference parameters reveal the formation of preferen- 
tial bonding within the same sublattice pairs. This phase 
can be viewed as a temperature evolution of the zero tem- 
perature LD structure, characterized by p^D = Cld = ^ 
and Pld = Cld ~ ^- Thermal fluctuations introduce 
defects in the diamond H bond network on the AB sub- 
lattices, and begin to populate the CD sublattices, which 
are rigorously empty at zero temperature. This phase is 
fourfold degenerate, and the alternative solutions can be 
obtained from Eq. (|20p by performing all possible circular 
permutations. 

In certain conditions, we also observe a peculiar phase, 
which we denote as "modulated" fluid (MF) phase, char- 
acterized by the following parameters 



p^ p^ p^ p^' 

e e-" e 



p'mf p'mf Pmf Pmf 















(21) 



Assuming p{^p > p^p, the two sublattices A and C are 
more populated than B and D (whence the term "modu- 
lated" ) , while there is no preference between the two arm 
configurations (whence the term "fluid" ) . This phase can 
be viewed as a temperature evolution of a zero tempera- 
ture state with p'-^-p = 1 (two fully populated sublattices) 
and pj^p — (two empty sublattices), or vice versa. Note 
that the two populated sublattices are next nearest neigh- 
bors, so that there is no interaction energy in these states. 
One can verify that these can actually be ground states 
of the GBHB and BL models, only if the two-body re- 
pulsive interaction is larger than the H bond energy, i.e., 
e/ri < -1. 




(b) 



0.8 1.0 



T/r] 



FIG. 5: GBHB model phase diagram, e/rj = —0.5: pres- 
sure vs temperature (a); density vs temperature (b). F, 
LD, HD denote the corresponding phases (see the text); dou- 
ble labels denote coexistence regions (b). Solid lines denote 
first-order transitions (a) or boundaries of coexistence regions 
(b); dashed lines denote second-order transitions; (thin) dash- 
dotted lines denote the TMD locus. Symbols display Monte 
Carlo data from Ref . [25l: solid squares denote first-order tran- 
sitions (a) or coexistence boundaries (b); open squares denote 
the TMD locus. 



Let us finally remark that, except the fully homoge- 
neous F phase, all the phases described above (HD, LD, 
and MF) are symmetry-broken phases, which do not pos- 
sess the full translational symmetry of the bcc lattice. 



B. GBHB model 

In Fig. O we report the phase diagram of the GBHB 
model, for e/i] = —0.5. This ratio corresponds to the pa- 
rameters chosen by Girardi and coworkers in the original 
paper [2^. In the T-P diagram (Fig. [5^), we observe 
three different first-order transition lines. The first one 
separates the F and LD phases, whereas the other two 
occur between the LD and HD phases. All these lines 
are mapped onto coexistence regions in the T-p diagram 
(Fig. O)). Both the LD-HD first-order transition lines 
terminate in tricritical points, which are connected by a 
second-order line (i.e., a line of critical points) enclos- 
ing the LD phase. Another critical line separates the F 
phase from the HD phase, and terminates in a critical 
end-point. The LD phase exhibits a density anomaly, 
namely, a temperature of maximum (or minimum) den- 
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sity (TMD), at constant pressure. In the T-P diagram, 
the TMD locus is bounded between a minimum and a 
maximum pressure. The portion of the TMD hue joining 
these two points corresponds to density maxima, whereas 
the remaining two branches of the hne to density minima. 

Fig. [5] also reports some data points obtained by the 
Monte Carlo study of Ref. [25]. We find a remarkably 
good agreement both for first-order transitions (transi- 
tion lines and coexistence regions in the T-P and T-p di- 
agrams, respectively), and for the TMD locus. Such an 
agreement suggests that the tetrahedral cluster approx- 
imation is able to take into account the most relevant 
correlations present in the system. Let us recall that in 
Ref. [2^ the authors could not detect second-order tran- 
sitions, and interpreted the tricritical points as critical 
points terminating first-order transitions between homo- 
geneous fluid phases of different densities. The misin- 
terpretation was probably due to the fact that no order 
parameter was investigated. Indeed, in Ref. 27 we col- 
lected several evidences supporting the existence of the 
ordered phases predicted by the cluster approximation, 
and the validity of the corresponding phase diagram. 

The previous results show that unfortunately the phase 
diagram of the GBHB model is quite far from that of real 
water. The main difficulty is that the density anomaly 
appears in a region where the system still exhibits orien- 
tational order (i.e., the icelike LD phase). Furthermore, 
the LD phase is never less dense than the corresponding 
fluid phase, as shown by the slope of the F-LD transition 
line in the T-P diagram, and no liquid-vapor coexistence 
is observed. We have explored the possibility of removing 
at least the first difficulty, by reducing the relative im- 
portance of the attractive H-bond interaction (which is, 
the one responsible for orientational order) , with respect 
to the repulsive close-packing interaction. 

In Fig. [HI we report the phase diagram for e/77 = —0.8. 
It turns out that the stability of the LD phase with re- 
spect to temperature is actually reduced, but our diffi- 
culty is not solved, since the TMD locus is displaced to- 
ward lower temperatures, remaining inside the LD phase 
region. Furthermore, a modulated fluid (MF) phase 
becomes stable at finite temperature, which makes the 
phase diagram even more complex. The MF phase is 
separated from the ordinary fluid (F) phase by a second- 
order transition and from the LD and HD phases by two 
different flrst-order transitions. The physical mechanism 
underlying the onset of the MF phase can be qualitatively 
understood by energetic arguments. On the one hand, 
the increased importance of nearest-neighbor repulsion 
favors conflgurations in which occupied sites have empty 
neighbors (by the way, this also enhances stability of the 
LD phase with respect to pressure). On the other hand, 
the reduced importance of H bonding favors thermody- 
namic states with no preference for particular molecule 
orientations. 

We have also investigated the opposite regime, in which 
the H-bond interaction dominates with respect to the 
close-packing interaction (i.e., 77 ^ |e|). We have ob- 
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FIG. 6: GBHB model phase diagram, e/r; = —0.8: pressure 
vs temperature (a); density vs temperature (b). Labels and 
lines are defined as in Fig. [S] MF denotes the modulated fiuid 
phase. 



served that the pressure and density ranges on which the 
LD phase is stable become smaller and smaller, but the 
phase diagram does not exhibit signiflcant changes with 
respect to the case tji] = —0.5. 



C. BL model 

As mentioned in the Introduction, the BL model can be 
viewed as a version of the GBHB model, with asymmetric 
bonds. The interaction parameters chosen by BL in the 
original paper [l8l | , based on a flt to real thermodynamic 
properties, correspond to e/rj = —0.08, which is, to a 
regime of dominating H-bond interaction. The phase di- 
agram we have obtained for this case is reported in Fig.[7l 
This phase diagram is, both qualitatively and quantita- 
tively, very similar to the one of the GBHB model for 
the same e/rj value. We just observe a sort of "rescaling" 
of the whole phase diagram toward lower temperatures. 
This effect is of entropic nature, and is related to the 
presence of the sign conflgurations. 

Apart from these changes, the most remarkable fact 
in the above phase diagram is that we flnd no trace of 
the vapor- liquid coexistence, claimed by BL [l^. We 
have already mentioned that the reason of such a dis- 
crepancy is a homogeneity hypothesis, on which the cal- 
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FIG. 7: BL model phase diagram, e/rj = —0.08: pressure 
vs temperature (a); density vs temperature (b). Labels and 
lines are defined as in Fig. [5] 



culations of the cited paper are based. As shown in the 
previous Section, our approach aUows one to impose an 
analogous homogeneity constraint in a straightforward 
way. We have repeated the phase diagram calculation, 
taking into account this constraint, and the results are 
reported in Fig. [51 Of course, we no longer observe the 
ordered phases, but only homogeneous (fluid) phases. At 
low temperature, there appear two first-order transition 
lines (and related coexistence regions), allowing us to dis- 
tinguish three phases with different densities. In particu- 
lar, the transition between the lowest density phase and 
the intermediate density one might be identified with the 
vapor-liquid transition. The other transition, between 
the intermediate and the highest density phases, reminds 
us of the conjecture about the liquid-liquid transition in 
metastable water, put forward by Stanley and cowork- 
ers [13]. Unfortunately, a more detailed analysis shows 
that this "homogeneous" phase diagram carries several 
inconsistencies. First of all, one can observe that the 
entropy is negative below a certain temperature, which 
depends on pressure. The TMD locus and most part of 
the supposed liquid-liquid coexistence region (including 
the critical point), lie below the zero-entropy line. More- 
over, making use of Eq. pip , it is possible to compute 
the continuation of the F-HD critical line in the stability 
region of the LD phase. This line marks a boundary, be- 
yond which the free-energy minimum associated to the 
homogeneous solution becomes a saddle point, i.e., the 
homogeneous phase becomes thermodynamically unsta- 
ble. As shown in Fig. [51 such a condition comes true in 



FIG. 8: BL model phase diagram in the homogeneity as- 
sumption, e/rj = —0.08: pressure vs temperature (a); density 
vs temperature (b). Dashed and dotted lines denote respec- 
tively the stability limit and the zero-entropy locus. Other 
lines are defined as in Fig. [5l 



a very large portion of the phase diagram. By the way, 
let us also remark that, although for brevity we do not 
report the results, totally analogous phase diagrams have 
been obtained for the "homogeneous" GBHB model. 

We have also reproduced BL's original calculation p^ . 
in order to compare the results with those obtained by 
our approach, in the homogeneity assumption. Fig. [9] 
reports the vapor-liquid binodals and some coexistence 
lines in the density-temperature and density-entropy 
planes. It turns out that, for not too high densities, the 
two methods yield very similar results. Conversely, a 
relevant discrepancy appear in the high density regime, 
which is, the phase diagram computed by BL [3| does 
not exhibit the "liquid-liquid" transition at all. Such 
a discrepancy is ultimately due to the different clus- 
ter considered for the approximation, namely, a nearest- 
neighbor pair in BL's calculation and the tetrahedral 
cluster in the current one. The (pseudo) liquid-liquid 
transition is indeed a reminiscence of the distinction be- 
tween the LD and HD structures, which the tetrahedron 
approximation is able to account for (although, in the 
homogeneous assumption, only at the level of local cor- 
relations), at odd with the pair approximation. 

Let us finally remark that, according to Fig. [SJd, even 
BL's solution turns out to suffer from the negative- 
entropy problem. More generally, one can see that all 
the aforementioned inconsistencies do not depend on the 
type of approximation (except the fact that the pair ap- 
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FIG. 9: Liquid-vapor coexistence region for the BL model in 
the homogeneity assumption, e/r; = —0.08: temperature vs 
density (a); entropy per molecule vs density (b). Thicker and 
thinner lines refer to our tetrahedral cluster approximation 
and BL's pair approximation, respectively. Solid lines denote 
the binodals; straight dotted lines denote coexistence lines at 
different temperatures. 



proximation does not reveal the liquid-liquid transition) 
but rather on the artificial homogeneity constraint. 



D. Bell model 

The Bell model is quite different from the GBHB and 
BL models, as the nearest-neighbor interaction is attrac- 
tive, while the close-packing repulsion effect is modeled 
by a three-body interaction (see Section II). For the pa- 
rameter set of the original paper [ll[ , we have obtained 
the phase diagram reported in Fig. [TOl We can indeed 
observe several differences, with respect to the phase di- 
agrams presented above. First of all, there appears a co- 
existence region between two homogeneous (fluid) phases 
at different densities, which we can identify with a vapor 
and a liquid phase, respectively. This new feature is likely 
to be related to the quite strong orientation-independent 
attractive interaction (e > 0) , which is absent in the other 
models. The ordered LD and HD phases are still present, 
but with clearly different phase transitions. The F phase 
turns out to be more stable with respect to pressure, so 
that a direct F-LD transition takes place, at odd with 



FIG. 10: Bell model phase diagram, ejr] — 2, S'y/rj — —5/4: 
pressure vs temperature (a); density vs temperature (b). La- 
bels and lines are defined as in Fig. [5l Different fluid phases 
with lower and higher densities are denoted by Fvap and Fuq, 
respectively. 



the GBHB and BL models. Such a transition is partially 
flrst-order (at lower pressure) and partially second-order 
(at higher pressure). As a consequence, a multicritical 
point appears, at which three critical lines (F-LD, LD- 
HD, and F-HD) merge. Density anomalies can still be 
observed in the LD phase. The TMD locus is made up 
of two different branches: a lower pressure one, corre- 
sponding to density maxima, and a higher pressure one, 
corresponding to density minima. The F phase does not 
display any density anomaly. 

Let us now consider the phase diagram obtained impos- 
ing the homogeneity constraint (Fig. [TT|) . Of course, no 
ordered phase is present, whereas the vapor- liquid coexis- 
tence region continues down to zero temperature. Let us 
note that Bell's original calculation [ll| is equivalent to 
the tetrahedral cluster approximation for homogeneous 
phases, so that the vapor-liquid binodals, obtained by 
the former, turn out to be exactly superimposed to those 
reported in Fig. 1111 In the very low temperature re- 
gion, we also observe a "liquid-liquid" coexistence, like 
in the BL and GBHB models. The latter feature was 



not pointed out by Bell [llj , most likely because at that 
time the exploration of supercooled water physics was at 
the very beginning [ij], and the "second critical point" 
conjecture had not been proposed yet [3^ . Anyway, the 



11 



1.5 ■ 



(a) 




T/ri 



FIG. 11: Bell model phase diagram in the homogeneity as- 
sumption, e/rj = 2, 37/77 = —5/4: pressure vs temperature 
(a); density vs temperature (b). Lines are defined as in Fig.[8l 



entropy computed by the cluster approximation turns out 
to be negative in this region. Conversely, the TMD locus 
lies mostly in the positive entropy region, and exhibits 
the correct experimental slope. The latter was indeed 
one of the most striking results of the Bell model. Unfor- 
tunately, the stability limit, which in this case is made 
up of both the F-HD and F-LD second-order transition 
lines (and a metastable continuation of the latter), sug- 
gests that, even for the Bell model, all the interesting 
anomalies take place in a region where the homogeneous 
F phase is thermodynamically unstable. 

We have compared our results for the Bell model with 
some Monte Carlo data obtained by Whitehouse and 
coworkers [l3|. In Fig. [T^] we report the density as a 
function of the chemical potential at fixed temperature 
T/t] = 0.5. At a first glance, it turns out that the tetrahe- 
dral cluster approximation does not fit the Monte Carlo 
results so closely as for the GBHB model. This fact is 
probably to be ascribed to stronger correlations, present 
in the Bell model, arising from the three-site interaction. 
We expect that the observed discrepancy is not so rele- 
vant to affect the qualitative structure of the phase dia- 
gram. Anyway, it is evident that the symmetry-broken 
solution (solid line) is much closer to the Monte Carlo 
results than the homogeneous one (dashed line). In par- 
ticular, the LD-HD symmetry change (see Subsection A) 
provides an explanation for the transition observed at 
higher chemical potential. Therefore, we argue that the 
simulations of Whitehouse and coworkers [ij] actually 
describe ordered phases of the same type predicted by 
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FIG. 12: Density vs chemical potential at constant temper- 
ature T/77 = 0.5 for the Bell model, e/r? = 2, 87/77 = -5/4. 
The dashed and solid lines denote the predictions of the 
cluster-variational calculation, respectively with or without 
the homogeneity constraint. Symbols denote Monte Carlo re- 
sults from Ref. |l3| (the thin dotted line is an eye-guide) . 



our approximation. This fact was not recognized in the 
cited paper, most likely because the authors did not in- 
vestigate any kind of order parameter. 



VI. CONCLUSIONS 

The basic ideas of the present paper originated in a 
previous one [2^ . in which we showed that a gener- 
alized first-order approximation on a tetrahedral clus- 
ter was able to reproduce several Monte Carlo results 
for a waterlike network-forming lattice model (GBHB 
model) 25]. The approximate theory also predicted that 
the two denser phases, identified as liquid in Ref. [25l . 
were in fact characterized by two different ordered H- 
bond-network structures. In the current work, we have 
exploited the same cluster approximation, to extend the 
analysis of the GBHB model to different parameter values 
and to revisit other models with similar features (tetra- 
hedral molecules, bcc lattice), already known in the lit- 
erature (BL and Bell models). 

We find that analogous ordered phases are present in 
all the models investigated. As mentioned in the text, 
several previous investigations did not point out these 
phases. In particular, while Girardi and coworkers did 
not re cog nize order-disorder transitions in their simula- 
tions [25|, as they did not compute order parameters, 
other studies ignored the ordered phases by imposing ho- 
mogeneity in an analytical way [ll|, [H, [l3| . Our results 
show that the onset of orientational order produces sub- 
stantial changes in the phase diagrams, which turn out to 
be very complex and rich, but quite far from real water. 

Let us recall that, in principle, the ordered phases are 
an interesting feature of the current models, since they 
might be regarded as a description of two different ice 
forms. Unfortunately, such an interpretation gives rise 
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to some difBculties. The most relevant one is that the 
icehke LD phase exhibits a temperature of maximum den- 
sity, which would be indeed expected in the liquidlike F 
phase. Conversely, the F phase turns out to be rather 
trivial, as it does not exhibit any density anomaly, and, 
in the GBHB and BL models, not even a vapor-liquid 
transition. These difficulties seem to be somehow related 
to the fact that, in these models, orientational order is 
exceedingly stable with respect to thermal fluctuations. 

We have explored the possibility that a suitable choice 
of the parameter set could yield a more realistic phase 
diagram, in the framework of the same models. Sev- 
eral failed attempts have led us to believe that this is 
likely not to be the case. Indeed, the physical mecha- 
nism underlying density anomalies (which is, correlation 
between low energy and high specific-volume configura- 
tions) is also responsible for enhancing thermodynamic 
stability of the ordered LD phase. Accordingly, a param- 
eter change oriented, for instance, to lowering the LD-F 
transition temperature, usually turns out to lower the 
temperature of maximum density as well. In Section |VB] 
we have reported an example of such an effect for the 
GBHB model. 

Let us remark that analogous arguments might prob- 
ably hold for other similar models, not studied in this 
paper. For instance. Bell and Salt [13] have shown 
that their model predicts open- and close-packed ice- 
like phases, which are respectively equivalent to the LD 
and HD phases of the current models. More recently, it 
was noticed ^l| that the liquidlike phase of the Roberts- 
Debenedetti model is metastable at "ordinary" temper- 
atures. Indeed, the authors of Ref. ^2^} devise particular 
techniques to prevent the Monte Carlo dynamics from 
"falling" into ordered states. Two authors of the present 
paper have verified the onset of analogous ordered states 
in a simplified (symmetric bonds) version ^23] of the 
Roberts-Debenedetti model, although this result has not 
been published yet. All these facts suggest that a quite 
large class of waterlike lattice models should be critically 
reconsidered. 

Concerning the stability of the H-bond network struc- 
tures, let us note that, in the models under investigation, 
this is necessarily overemphasized by the regular bcc lat- 
tice, which imposes a strong directional correlation on 
the bonds. Such a correlation is not present in the real 
(off-lattice) system. In particular, it is known from ex- 
periments that directional correlation in real water is al- 
most completely lost after the second consecutive bond 
(see Ref. [l] and references therein). Therefore, a more 
realistic model might be provided by a suitable random 
lattice, preserving the local geometric structure of the 
H-bond network (up to the second consecutive bond). 
We expect that such a lattice might hinder the onset of 
orientational order, though retaining the essential physi- 
cal properties of the liquidlike phase. Of course, even a 
random-lattice description involves a simplification, be- 
cause, in principle, directional disorder should not be an 
a-priori assumption, but rather a prediction of the theory. 



In any case, a very simple realization of the random 
lattice, satisfying the requirement on the local tetrahe- 
dral geometry, is the Husimi lattice made up of tetra- 
hedral blocks, which we have mentioned in Section IV 
(Fig. |4]). As noted there, the tetrahedral cluster approx- 
imation happens to be exact for this particular lattice. 
Furthermore, it is known that, in general, Husimi and 
Bethe lattices are characterized by a random graph struc- 
ture, which is locally treelike, but contains closed loops 
of different (even and odd) lengths. The latter fact can 
produce frustration, which is expected to forbid periodi- 
cally ordered states [i^ like those appearing on a regular 
lattice. As a consequence, it turns out that the different 
models treated in this paper, properly redefined on the 
Husimi lattice, cannot exhibit ordered phases at all. The 
exact thermodynamic equilibrium states for these mod- 
els are given by the "homogeneous" stationary points of 
the appropriate variational free energy. This argument 
allows one to reinterpret the "homogeneous" phase dia- 
grams, presented in this paper, which exhibit interesting 
similarities with real water thermodynamics. In the new 
interpretation, such results are no longer artifacts of the 
homogeneity assumption, but rather exact solutions for 
the corresponding Husimi lattice models. 

At first sight, the above idea might look a bit tricky, 
but in fact a random lattice is not much more arbitrary 
than any regular lattice, for modeling a fiuid. Both types 
of lattice should be regarded as simplifying assumptions, 
which affect the model predictions in opposite ways. In 
the regular lattice case, the stability of the ordered phases 
is overestimated, whereas, in a generic random lattice, 
it is underestimated. In the simplest case of a Husimi 
lattice, whose exact solution is almost analytically avail- 
able, periodically ordered states turn out to be rigorously 
forbidden. In other words, the price paid for simplic- 
ity is that the model is no longer able to describe ice- 
like phases. Let us also remind that a Husimi lattice is 
an infinite-dimensional system, so that all its properties 
have necessarily a mean-field nature. As a consequence, 
the model cannot yield realistic estimates of critical ex- 
ponents, which turn out to belong to the mean-field uni- 
versality class. 

Let us also stress the fact that the tetrahedral cluster is 
not the only possible choice for the building block of the 
Husimi lattice. This choice turns out to be quite effective 
for the current models, as it is able to distinguish two dif- 
ferent (low- and high-density) local molecular packings. 
As mentioned in Section IV C[ these two packings are not 
only the elementary structures of the ideal LD and HD 
networks (which are not relevant on the random lattice), 
but they are also expected to play a role in the possible 
onset of a liquid-liquid phase transition. 

Let us finally note that, in the light of the new in- 
terpretation, the negative-entropy phenomenon suggests 
that the waterlike Husimi lattice models might predict a 
replica symmetry breaking, in analogy with the models 
studied by Mezard and coworkers [i^, |4^. This possi- 
bility might even be relevant for describing the glassy 
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states of water (amorphous ices), which have been stud- 
ied intensively in the last years by both experiments and 



simulations [41] . We are going to investigate this issue in 
a future work. 
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